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ABSTRACT 

An  absolute  scalar  pseudo-viscosity  ("q")  expression  for  handling  shocks 
automatically  in  numerical  calculations  of  three-dimensional  fluid  dynamics 
problems  is  derived.  The  expression  is  a  generalization  of  that  derived  by 
von  Neumann  and  Richtmyer  in  the  sense  that  the  expressions  are  identical 
for  a  plane  shock  wave.  The  basic  assumption  made  is  that  a  shock  will  be 
formed  in  a  material  when  there  is  compression  in  the  direction  of  accelera¬ 
tion;  this  assumption  seems  "reasonable"  physically,  since  this  condition 
indicates  that  a  grid  established  about  a  point  with  such  compression  will 
collapse,  i.e.,  the  "rear"  of  the  grid  will  overtake  the  "front.  " 

The  expression  "q"  is  derived  for  a  generalized  curvilinear  coordinate 
space  and  the  expression  is  given  for  the  three  commonly  used  coordinate 
systems.  The  expression  is  also  written  for  the  plane  rectangular  case  and 
the  cylindrical  case  with  no  angular  motion  to  be  used  in  connection  with  the 
HEMP  code. 


INTRODUCTION 

The  use  of  a  pseudo-viscosity  to  handle  shocks  automatically  in  numer¬ 
ical  calculations  of  fluid  dynamics  problems  is  well  known.1  The  first  programs 
written  were  lor  the  solutions  of  problems  in  one  space  dimension  for  which 
the  pseudo- viscosity  ("q")  devised  by  von  Neumann  and  Richtmyer  was  derived. 
However,  the  expression  for  "q"  in  this  one -dimensional  geometry  did  not 
carry  over  directly  into  more  than  one  dimension.  One  of  the  expressions 
_ 

von  Neumann  and  R.  D.  Richtmyer,  "Finite  Difference  Methods  for  Initial- 
value  Problems"  (Interscience  Publishers  Inc.,  New  York,  1957)  p.  205  ff. 
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used  for  a  two  dimensional  program  involves  the  time  rate  of  change  of  rela- 

§  O  #  f 

tive  volume  divided  by  relative  volume  (V/V)  which  in  one  dimension  is 
equivalent  to  the  gradient  of  the  velocity  used  by  von  Neumann  and  Richtmyer. 
However,  in  a  cylindrical  or  spherical  implosion,  the  term  V/V  will  become 
negative  (forming  a  shock)  simply  because  the  mass  is  being  forced  into  the 
center  even  though  the  rear  of  a  particular  grid  is  not  overtaking  the  front  of 
that  grid.  The  derivation  given  here  overcomes  this  difficulty  and  is,  in 
fact,  a  straight  forward  generalization  of  the  expression  devised  by  von 
Neumann  and  Richtmyer  in  the  sense  that  the  "q"  was  required  to  indicate 
that  compression  is  present  in  the  direction  of  acceleration,  to  turn  itself  on 
when  this  condition  develops,  and  to  give  the  same  expression  when  a  plane 
shock  passes  through  a  Mone  dimensional  material.'1  The  expression  derived 
differs  from  theirs  in  that  it  is  formed  from  vector  quantities  in  three  space 
dimensions  and  the  plane  shock  is  not  required  to  be  parallel  to  one  of  the 
coordinate  planes  in  order  for  the  nq,!  to  work.  The  basic  assumption  made 
is  that  a  shock  will  be  formed  when  there  is  compression  in  the  direction  of 
the  acceleration.  The  Mq"  so  derived  is  an  absolute  scalar  formed  from  the 
product  of  an  absolute  second-order  tensor  and  two  vectors  and  could  be  used 
even  in  three  dimensions. 

THE  DERIVATION  OF  A  GENERALIZED  VON  NEUMANN  SCALAR  "q" 
WITH  DIRECTIONAL  PROPERTIES3 

Consider  a  body  that,  under  some  unspecified  force  system,  is  under¬ 
going  motion  which  will,  in  general,  include  rigid  body  translation  and 
rotation  as  well  as  deformation  (Fig.  1).  Let: 

XK(K  =  1,2,3)  refer  to  the  coordinates  of  a  point  in  the  body  in 
a  fixed,  Lagrange,  generalized  coordinate  system. 

xk(k  =  1,  2,  3)  refer  to  the  coordinates  of  the  same  point  in  a 
moving,  Eulerian,  generalized  coordinate  system. 

Z^(K  =  1, 2,  3)  be  a  fixed  rectangular  coordinate  system. 
zk(k  =  1,  2,  3)  be  a  moving  rectangular  coordinate  system. 

^Mark  L.  Wilkins,  "Calculation  of  Elastic-Plastic  Flow,"  Lawrence 
Radiation  Laboratory  (Livermore)  UCRL-7322  (1963). 

o 

°The  derivation  and  notation  are  taken  largely  from  Chapters  1  and  2  of 
Nonlinear  Theory  of  Continuous  Media  by  A.  Cemal  Eringen  (McGraw-Hill 
Book  Co.,  Inc.,  New  York,  1962). 
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12  3  K 

P  be  the  position  vector  to  (X  ,  X  ,X  )  in  the  Z  system. 

12  3k 

p  be  the  position  vector  to  (x  ,  x  ,  x  )  in  the  z  system. 

b  be  the  vector  that  represents  the  rigid  translation  of  the  body. 

'v  K  1  2  3 

G. ,  be  a  base  vector  in  the  direction  of  X  from  (X  ,  X  ,  X  ). 

—  L  k  1  2  3 

g.  be  a  base  vector  in  the  direction  of  x  from  (x  ,  x  ,  x  ). 

K 

t  be  the  time. 


Fig.  1 

Then  the  motion  at  any  time  is  expressed  as: 

VK  ,rK,  12  3.,  k  k/vl  „2  v3  ..  /lX 

X  =  X  (x  ,  x  ,  x  ,  t)  or  x  =  x  (X  ,  X  ,  X  ,  t).  (1) 

From  Eq.  (1)  we  can  determine  the  path  taken  by  a  point  in  the  unde- 
-1-2-3  k 

formed  body  at  (X  ,  X  ,  X  )  as  a  function  of  xK  and  t  and  conversely,  we  can 

12  3  K 

determine  the  path  taken  to  (x  ,  x  ,  x  )  as  a  function  of  X  and  t  (in  this 
treatment  the  inverse  of  X^(x\  x2,  x3,  t)  is  assumed  to  exist).  In  other  words, 
given  a  point  in  the  undeformed  body,  we  can  determine  where  it  goes  and 
given  a  point  in  the  deformed  body  we  can  determine  where  it  came  from  in 
the  undeformed  body. 

Consider  a  differential  change  in  P  in  the  generalized  Lagrange  coordinate 
system.  This  differentia]  is  expressed  by: 

dP  =  G„  dX1^  =  G.  dX1  +  G„  dX2  +  G„  dX3  (2) 


and  the  square  of  the  magnitude  of  this  vector  is: 


where  Gj^ 


dS2  =  dP  •  dP  =  G^.  dXK  dXL, 

~  ~  KL 

-  Gt^  •  Gt  is  the  metric  tensor  for  the 
~,K  ^  L 


(3) 

system  of  coordinates. 
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Similarly,  for  the  Eulerian  coordinate  system  one  obtains: 


dp  =  g,  dx 

*%/ 


(4) 


along  with  its  magnitude  squared: 


where  g^ 


2  k  £ 

ds  =  dp  •  dp  =  gM  dx  dx  , 

k 

=  g  •  g  is  the  metric  tensor  for  the  x  coordinate  system. 


(5) 


The  quantities  dS2  and  ds2  can  be  interpreted  as  the  square  of  an 

element  of  arc  length  before  deformation  and  after  deformation  respectively. 

It  is,  therefore,  possible  to  obtain  a  measure  of  the  deformation  produced 

2  2 

by  the  motion  by  calculating  ds  -  dS  . 

The  material  or  substantial  derivative  is  a  time  derivative  with  the 
Lagrange  coordinate  XKheld  constant.  It  is  denoted  by  D/Dt  and  it  can  be 
shown  that  (see  Appendix  A): 


JD 

Dt 


(ds2) 


2  d,  „  dx^  dx^ 
k  SL 


(6) 


where  d  ~  \  (v.  .  v.  ,  )  is  the  symmetric  part  of  the  velocity  gradient,  is 

iv£  k;X 

called  the  deformation  rate  tensor,  and  is  an  absolute  second-order  tensor. 


v  „  is  the  covariant  partial  derivative  of  v,  (see  Appendix  A,  Eq.  (A.  1)). 
k;i  K 

From  (6)  we  see  that: 


1  _D 
2  Dt 


ds 


(ds2)  =  ds  Dt  (ds)  =  2dk SL 


,  k  ,  i 
dx  dx 


ds  ds 


or 


1  D  ,  ,  , 
ds  Dt  (ds) 


,  k  JL 
dki  n  n 


=  d 


(n) 


(7) 


where  dx^ /  ds  =  n^  is  the  unit  vector  along  ds,  and  is  called  the  rate  of 
stretching  or  rate  of  deformation. 

Since  dS  is  a  function  of  X  only: 


^  (ds  -  dS)  =  (ds)  (8) 

so  we  see  that  (7)  is  a  measure  of  the  rate  of  change  of  deformation  with 

respect  to  time  divided  by  the  current  length  of  an  element  in  the  direction 

n.  Hence,  (7)  can  be  used  to  determine  whether  there  is  tension  or  com- 

pression  in  the  direction  n,  where  n  is  an  arbitrary  unit  vector.  In  other 

~  .12  3. 

words,  given  an  arbitrary  direction  and  a  point  p(x  ,  x  ,  x  )  in  a  deforming 
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body,  (7)  can  be  used  to  determine  whether  an  element  of  length  in  that  direc¬ 
tion  is  lengthening  or  shortening  with  time  at  a  selected  time  t. 

Consider  a  body  in  motion,  select  a  point  in  this  body  at  p(x  ,  x  ,  x  ). 

One  of  the  characteristics  of  a  von  Neumann  "q"  is  that  it  is  activated  if  the 
motion  of  p(x\x^,x^)  is  a  deceleration,  i,e.,  the  grid  is  collapsing  in  the 
direction  of  acceleration.  It  is  clear  that  if  the  vector  n  in  (7)  is  a  unit  vector 
in  the  direction  of  the  acceleration  vector,  then  (7)  will  be  able  to  measure 
this  effect,  i.e.  ,  if 


d(n)  =  ds  Dt  (ds)  <  °* 

then,  with  nK  -  A  /|A|,  the  grid  is  collapsing  in  the  direction  of  motion, 

where  A  is  the  acceleration  vector. 

It  is  important  to  note  here  that  is  completely  independent  of  the 

coordinate  system  that  is  used  to  express  d^  and  n  .  This  is  easily  seen 

from  the  method  of  derivation  since  the  element  ds^  =  dp  •  dp  is  a  scalar 

derived  from  the  scalar  product  of  two  vectors.  This  can  also  be  seen  from 

(7)  as  follows;  Since  d.  is  an  absolute  covariant  tensor  of  second  order,  a 

k*  k  —  k 

transformation  of  coordinates  from  x  to  x  through 


dx 


3x 

a  1 
o  x 


dx 


results  in  a  transformation  of  d^  to  dj^  by 


d! 


Q  m  n 
ox  ox 


kjf 


inn  g— k  g—  & 


k.  —  k 

by  definition.  The  contravari ant  vector  n  will  transform  to  n  by 


n 


m  3x 


k 


n 


9x 


m 


We  now  form,  using  these. 


dh 

1 

1 

8xm 

9xn 

n  n  =  ci 

mn 

9xk 

ax'2 

db 

—  k  —  jg 

8xm 

8xk 

n  n  =  d 

mn 

axk 

3xr 

d|  „ 

— k— je 

e  m 

c  n  r 

n  n  =  d 

6 

6  n 

ki 

mn 

r 

P 

r  3x 


n 


3x 


9xn  ax-2 


ax'2  9xp 


p  dx* 

axp 

n~  nP 
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where 


1  if  i_=  j 


0  otherwise 


=  Kronecker  delta 


or 


— k  —St  ,  m  n 
d'  n  n  =  d  n  n  . 
ki  mn 


This  invariance  with  respect  to  an  arbitrary  coordinate  system  is  the  chief 
property  of  interest  in  this  expression.  In  the  HEMP  code  the  "Q"  is  based 
on  V/V,  where  V  is  the  relative  volume.  In  a  spherical  or  cylindrical  im¬ 
plosion  the  expression  V/V  becomes  negative  simply  because  the  mass  is 
being  crowded  in  to  a  smaller  and  smaller  space  while  in  the  direction  of 
acceleration,  the  front  surface  of  the  grid  (the  one  toward  the  center)  is 
moving  faster  than  the  rear  surface.  The  quantity  d^nj  would  indicate  this 
clearly,  if  £  were  a  vector  in  the  direction  of  the  acceleration  vector. 

If  we  consider  a  shock  moving  through  a  material  at  rest  and  making 

12  3 

an  angle  with  the  coordinate  system  x  ,  i.e. ,  not  along  any  of  x  ,x  ,x  ,  then 
a  coordinate  system  can  be  selected  such  that  one  coordinate  a  lies  along  the 
acceleration  vector  and  the  other  two  coordinates  (/3,Y)  are  normal  to  a. 

Then 

n*  =  1  ,  n^  =  0  ,  n^  =  0 

and  ,  9v 

d/  \  -  d/  x  -  d  +  v  )  =  ■  a  —  , 

(n)  (a)  aa  2  cr,a  a\a  da 

where  v  is  the  velocity  in  the  a  direction  and,  in  this  case,  is  in  the  direction 
a 

of  the  shock.  In  other  words,  d^  =  d^  is  the  partial  derivative  of  the  velocity 

in  the  direction  of  the  shock  witlTrespeet  to  the  space  coordinate  in  the  direction 

of  the  shock.  We  see,  therefore,  that  this  is  the  same  expression  used  by  von 

Neumann  and  Richtmyer'*'  in  the  "activating"  part  of  the  artificial  viscosity  and 

all  that  has  been  done  here  is  to  express  it  in  a  generalized  form. 

The  sticky  part,  of  the  problem  involves  the  selection  of  the  characteristic 

length  to  be  used  to  multiply  the  activation  term  and  spread  the  shock  over  the 

same  number  of  zones  in  the  numerical  solution,  von  Neumann  and  Richtmyer 

used  AX  the  Lagrange  grid  size  in  the  direction  of  the  shock  (see  Ref.  1, 

9 

p.  215  ff).  In  this  case  our  grid  is  established  in  the  system  X  ,  is  allowed  to 
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deform  in  the  system  and  is  not  in  general  in  the  direction  of  the  shock. 
Conceptually  a  grid  length  A  A  could  be  established  in  the  shock  direction, 
corresponding  to  a  above,  and  this  could  be  transformed  to  the  X  coordinate 
system,  but  in  practice  this  cannot  be  done  since  the  grids  are  already  estab- 
lished  in  the  X  system.  The  approach  taken  here  was  to  use  the  established 
grid  to  determine  a  representative  length  in  the  direction  of  the  acceleration. 

Two  lengths  were  tried  (AS1  and  ASg).  ASX  was  formed  as  the  arith¬ 
metic  average  of  the  projections  of  the  two  grid  diagonals  in  a  two-dimensional 
grid  on  the  acceleration  direction.  AS2  was  formed  as  the  square  root  of  the 
sum  (or  average)  of  the  geometric  averages  of  the  separate  components  of  the 

diagonals  projected  on  the  acceleration  direction.  If  the  figure  below  repre- 

1  2 

sents  the  grid  at  time  step  t  and  the  quantities  n  and  n  are  understood  to  be 
calculated  at  the  midpoint  of  the  grid  at  time  step  t,  then  the  expressions  used 
are: 


.4! 

1 _  J3 

r  ^ 

T 

2 

GRID 


AS^- 

|<X3  - 

as2  = 

[(a1)2 

|<X3  -Xjjn1  +-  (Y3  -  Y^n2!  +  |(X2  -  X4)  n1  +  (Y2  -  Y^n* 


2,2 


^1  1  1  2  4 


1  1  '  2  ~  4 


1/2 


where  |  B  |  is  the  absolute  value  or  magnitude  of  B,  and  where 


it/!1 


n 


1  =  P  *.) 


+  f1 


|i  t  i1  +  r 

'p 


17-  ‘'i‘  1'2 


are  the  general  expressions  for  the  direction  cosines  in  the  acceleration 
direction  (see  Ref.  3,  p.  104).  (Here  t1^  is  the  stress  tensor,  f1  the  body 
force  per  unit  mass,  and  p  the  mass  density.  )  The  factors  of  1/2  needed 
in  the  averages  are  to  be  embodied  in  a  multiplicative  constant. 

In  the  limit  of  a  shock  in  one  dimension,  let  Y  =  0,  then: 


ASX  = 

as2  = 


X 


3 


-XlHX2 
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So  the  grid  size  in  the  non-shock  direction  is  not  considered  in  the  one¬ 
dimensional  shock  as  would  be  desired.  Further,  if  the  grid  is  square,  then 
|X3  -  Xt  |  =  |X2  -  X4  |  and  ASj  =  AS2  =  AX  which  is  the  factor  used  by  von 
Neumann  and  Richtmyer. 

In  the  problems  worked  using  this  "q,"  there  was  no  obvious  difference 
between  the  two  AS's.  Some  experimentation  would  be  required  in  order  to 
make  a  selection. 

The  HEMP  code  is  a  two-dimensional  code  with  plane  or  cylindrically 
symmetric  geometry  and  the  expression  for  d^  is  (see  Appendix  B): 


ax 

ax 


(n1)2 


+ 


ay 

ay 


(  2x2 

(n  ) 


fax 

lay 


1 


n 


2 

n 


where  Y  is  either  the  radial  direction  or  the  other  rectangular  direction, 
this  code  we  have  used  both  a  linear  and  a  quadratic  "q"  as  follows: 
Linear: 


In 


q  =  0  if  d(n)  ^  0 

o  =  -  CL  T  d(n)  AS  ifd(n)<0 


Quadratic: 

q  =  °  if  >  0 
/v 

o  =  C0  4  (d(n)f  <AS)2  ifd(n)<0. 

where 

&nd  Cq  are  constants 
a  is  the  sound  speed 
Pp  is  reference  density 
V  is  relative  volume 
AS  =  ASj^  or  AS2  as  discussed  above 
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A 

x 


l  as 


XX 


3x 


+ 


9T  T 

_ *y+  .jar 

9y  y 


acceleration  in  x-direction 


p  =  mass  density 


acceleration  in  y-direction 


CONCLUSION 


The  expression  d^  derived  here  is  seen  to  be  a  measure  of  the  rate  of 
deformation  of  an  element  of  length  through  a  material  point  divided  by  that 
length  in  the  direction  n.  It  was  assumed  that  if  n  were  along  the  acceleration 
vector  through  the  point  and  if  d^  <  0  (compression),  causing  the  grid  estab¬ 
lished  about  that  point  to  collapse  so  that  the  rear  surface  is  overtaking  the 
front  surface,  then  a  shock  is  forming.  The  validity  of  the  "q"  depends  upon 
this  assumption  which  seems  to  be  the  assumption  made  by  von  Neumann  and 
Richtmyer  and  seems  to  be  "physically  reasonable."  d^  gives  an  exact  mea¬ 
sure  of  the  tension  or  compression  in  the  direction  n.  The  problem  is  to 
choose  the  "right"  direction.  The  direction  chosen  (acceleration  vector)  is  the 
one  which  indicates  that  particles  are  "bunching  up"  in  the  direction  of  impend¬ 
ing  motion,  which  is  assumed  to  indicate  the  formation  of  a  shock. 


These  terms  are  set  to  zero  for  a  plane  problem. 
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APPENDIX  A 


THE  MATERIAL  DERIVATIVE  OF  ds 

First  we  define  the  covariant  partial  derivative  of  a  vector,  which  is 

* 

indicated  by  the  semicolon  in  the  tensor  expressions  below  : 


r 


Ak 

k 

9x 

Im 

to 

9  A, 

V 

\  ;£ 

k 

=  a  a 

ox 

kJt 

L 

is  a  contravariant  v 

/  ^ 

k 

„2  n 
=  9  z 

9 

1m 

V,  1 

[  9xi  9xm  9 

A 


m 


A 


m 


(A.  1) 


n 


,k 


is  the  Christoffel  symbol  of  the  second  kind.  The  quantities  A  ^  and  A^.^,  are 
second-order  tensors.  Next,  we  define  the  material  or  substantial  derivative 
which  is  a  time  derivative  with  X^"  held  constant: 


.  „k  •  m  .  8B<  +  pk  .. 
-8T  +  Bi;mx 


m 


(A.  2) 


_  m  . 


with  Bk  =  B^x1,  x2,  x3,  t),  where  =  9xm/9t  -  Dxm/Dt  =  v“  is  the  velocity 
vector  in  the  Eulerian  system. 

12  3 

Now  the  quantity  g^  is  a  function  of  x  ,  x-  ,  x  only  and  by  Ricci's 


theorem: 


hence: 


gk  A; 


m 


=  0, 


D 


Dt  (gki)  ‘  °* 

By  considering  the  second  of  (A.  1): 

dxk  =  xk  „  dXK  at  any  time  t 
>  lx 

where  xkK  *  9xk/9XK  is  called  the  displacement  gradient. 


(A3) 


(A.4) 


"See  Ref.  3,  p.  43  9ff. 
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Now,  we  calculate  the  material  derivative  of  x  First  we  recall  that 

xkR  =  xk  K  (X1 ,  X2,  X3,  t)  but  by  the  first  of  (A.  1),  XK  =  X^x1 ,  x2,  x3,  t)  and 
hence  we  can  consider: 


xkK(X1,  X2,  X3,  t)  =  a^x1,  x2,  x3,  t) , 


(A.  5) 


then  using  (A. 2)  and  (A. 5): 


Dt  ^x,  K J  Dt  “  K 
Using  (A.  1)  in  (A.  6) 


3a 


k 


5-  =  D.  ak  (xl,x2,x3it)  =  ZJ<  +  ak, 


3ak  6  *  k 

D  k  K  u _ K  v  St  +  i 


Dt  x,  K  3t 


9x 


Ira 


at 


m  •  SL 
KX 


K;i 


k  . 


Now  the  ordinary  derivative  with  respect  to  time  of  a  ^is  given  by: 


d_ 

dt 


8i>kK+8akK.« 

~  a?  X 


•  jP  i 

and  is  as  noted  earlier  x  =  v  ,  hence  (A. 7)  becomes: 


D  k  d  I  k 


Dt  ,  K  dt.  I  K 


!  k  ) 


+ 


l  m  jP 
?  a  T,  v  . 


;j?m  ; 


K 


(A.  6) 


(A.  7) 


(A.  8) 
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If,  now,  we  think  of  a  K  =  x  ^(X  ,  X  ,  X  ,  t)  we  can  invert  the  order  of  differ- 

>  ^ 

entiation  in  the  first  term  on  the  right  of  (A. 8)  and  write: 


D  k  k  ,  I  m 

Dt  X,  K  V,  K  ‘  X,K 


(A.  9) 


k  k  1  2  3 

By  the  chain  rule  of  differentiation  considering  v  =  v  (x  ,  x  ,  x  ,  t)  and  using 
the  second  of  (A.  1): 


k  _  k  m 

,K  " 

so  that  using  (A.  1 ) 


V  =  V  X  jr 

,  K  ,  m  ,  K 


D  k 

x 


k  4- 
V  + 


Dt  ,  K  \  ,  m 


fk! 


,imj 


1  m  .  k  m 
v  x,  K  ‘  v  ;mx  ,K 


(A.  10) 


PC 

Since  D/Dt.  considers  X  constant,  multiplying  (A.  10)  by  dX  and  using  (A.4), 
yields: 

D  ,  k  k  i 
T5t  '  v  ;idx  • 


(A.11) 
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Now  we  are  in  a  position  to  consider  the  material  derivative  of  ds  . 
Using  (5): 


_D 

Dt 


,,  2,  D  f  ,  k  ,  Jt\  D  ,  ,  kv  ,  !L  .  A  k  D  ,  £. 

(ds  )  =  Dt  (gki  dx  dx  j  =  Dt  (dx  )  dx  gi \r0  dx  rff  <dx  ) 


Dt 


where  we  used  (A. 3)  to  remove  g^  from  inside  the  operator.  If  now  (A.  11) 
is  used: 


D  ,  ,  2.  k,m,i,  JL  ,  k  ,  m 

•prr  (ds  )  =  g.  „  v  dx  dx  +  g,  „  v  dx  dx  . 
Dt  &ki  ;m  &kf!  ;m 

❖ 


(A.  12) 


As  is  well  known  the  metric  g^  is  a  lowering  operator,  i.e.,  v^  -  g^v  and 
is  symmetric.  Using  this  and  changing  the  dummy  indices: 


Wt  (ds2)  '  (vi ;k  +  vk;J)dxkdx^  =  2  dW  dxkdx^ 


(A.  13) 


where  d^  r  77  (v^.^  +  v^.^)  symme^ric  Par‘t  °f  the  velocity  gradient  and 

is  called  the  deformation  rate  tensor. 


See  Ref.  3,  p. 


434. 
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APPENDIX  B 
EXPRESSIONS  FOR  d(n) 

A/ 

We  will  now  derive  the  expressions  for  d^  in  three  orthogonal  coordinate 
systems. 

Before  proceeding  it  is  necessary  to  introduce  the  concept  of  physical 
components  of  a  vector  and  a  tensor.  In  general,  the  components  of  a  vector 
or  a  tensor  do  not  have  the  same  units.  For  a  displacement  vector  with  con- 
travariant  components  referred  to  a  cylindrical  coordinates,  as  an  example, 
we  will  show  this. 

Let: 


u  r  u 


J?k 


be  the  displacement  vector.  The  base  vector  g^  is  defined  by: 


9p 


Mk  " 


3z 


m 


3x 


k 


3x 


p  ~m 


(B.  1 ) 


(B.  2) 


m 


In  a 


where  i  =  unit  vector  along  the  rectangular  coordinate  axis  z 

12  3 

cylindrical  coordinate  system  x  :  r,  x  -  0,  x°  =  z  in  the  usual  notation  and 
12  3 

z  -  x,  z  =  y,  z  =  z  in  the  rectangular  system  in  the  usual  notation  so  that: 


11  2  2 
z  =  x  cosx  ,  z 


1  .  5 

x  sin  x 


=  x 


from  which: 


(. 


dz 


m 


dx 


2  1.2 
cos  x  -  x  sin  x 


2 


sin  x 
0 


1  2 
X  cos  X 


\ 


1 


and 


2, 


gj  =  (cosx  )  i^  +  (sin  \  )  i2 
g0  ~  -(xhin  x2)i  +  (x1cosx2)ir 

W.L 

=  i3- 


Hence: 


=  1.  |g, 


=  x 


!  go  I  =  1. 


* 


See  Ref.  3,  p.  436  ff. 
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Now  if  u  is  a  displacement  vector,  it  has  the  units  of  length  L  and  since 
|g  |  -  1  and  |  |  =1,  u  and  uJ  have  units  of  length  L  but  |g2l  =  x  has  units 

L  whence  u  is  unitless.  Since  g^  =  g^  •  g^,  we  see  that 


(no  sum  on  k). 


(B.  3) 


In  order  to  express  the  components  of  a  vector  in  physical  components,  i.e., 
components  that  have  the  same  units  as  the  vector,  it  is  necessary  to  take 
parallel  projections  of  the  vector  on  unit  vectors  along  the  coordinate  curves. 
If  we  define  such  a  unit  vector: 


sk  =  Ik/jgkk  ' 


(k) 


then  we  define  the  physical  components  of  u,  denoted  by  u  by: 


(k) 

u  =  u  e 


k  ' 


(B.4) 


(B.  5) 


Since  a  vector  is  not  dependent  upon  the  coordinate  axis  used  to  measure  it, 
we  can  write: 

k  (k) 

U  Ik  =  U  £k 

and  using  (B.  4),  we  can  write: 

(k)  _  k  --  -  k  _  (k)/ 

u  -  u  N/% k  '  Ll  -  u  Ng kk  • 


(B.6) 


To  obtain  the  physical  components  of  the  covariant  components  we 
lower  the  index,  i.e. 


u.  =  g,  „  u 
k  6ki 


Ski  u 


U) 


Jt>  l  ' 


(B.  7 ) 


To  carry  this  definition  to  second-order  tensors,  a  tensor  B  is  converted  to 

k  a  ^ 

a  vector  B  by  the  use  of  a  vector  n  as  follows: 

-ok  ok 
B  =  B  ^n  , 

k  H 

then  using  (B.6)  on  B  and  n  one  obtains: 


B  (f)=B  W  gk  l  - 


(B.  8) 


It  will  be  understood  that  where  a  pair  of  repeated  indices  are  underlined, 
that  no  sum  is  to  be  taken. 
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Equation  (B.8)  gives  the  right  physical  components  of  B  The  left  physical 
components  are  given  by: 


B 


oo .  B  k 


Jekjg, 


(J>)  SL  ^  skk/6i  i 

In  general,  B(k^*  B(*k).  Now: 

,  .  k  &  ,  k  i 

d(„)  dMn  n  =  nkdi  n 

/V> 

using  (B. 6-B. 8): 


(B.  9) 


\n) 


(m) 


g 


f  j(k)  S  '  n(i)  j 


km 


m 


_i3_  z**. 

mjl  -J  ekk]\  Jg 


HU  J 


or 


d0->  =  II  Z,  Jgr 


®km 


mm  ®kk 


n(m)d(k)  nU) 
(£)  ’ 


k  H  m 

For  an  orthogonal  coordinate  system  g^m  =  0  for  k  1  m,  whence: 


d 


(n) 


n(k)d(k)  nU) 
n  d(i)  n 


(B.  10) 


and  writing  (B.10)  in  component  form: 

,(D  L (l) 


r  / 1  \  2  /  o  \  I  / n \  l  2  I o  \  r  / o  X  I  2 


Tn) 


d  (1)  •  *  +d 


:<2)  J"r<2) 

(2)Lr 


P)  L(3)l 
(3)  Ln 


n'  '  +  d  /0\  in’  , 


+  d(l()2)„<1)„(2>+d(2()1)n<1'n'2)+d(1>3)n(1>n(3) 
+  d(3>1)n(I>  „<:i>  >  d(2)3)n<2)  n(3)  +  d(3)2)n(2>  „<3). 


(B.  1 1 ) 


We  have  defined: 


dkjP  =  2  (vk;i  +  v*;k>  • 


and  we  obtain  from  this: 


,k  km  ,  1  (  k  ,  km 

di“8  dmi=2iv^+S  ve,mr 

where  gkm  is  the  metric  tensor  for  the  reciprocal  base  vectors  and  its 
matrix  is  given  by: 
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Now 


so  that 


z  km.  /  .-1 

g  )  ^§km 


n 

v„  =  g„  v 
ft  ;m  &jfn  ;m 


ft  2 


1  (  1 

2  (v; 


k  ,  n  mk 

+  p„  v  g 
&in  ;m & 


)• 


The  matrix  for  can  thus  be  obtained  from: 

(^H  f(V;kJ  +(8inv;m«mk)T)' 

T 

m  m 

where  (A  n)  indicates  the  transpose  of  the  matrix  (A  ). 

For  an  orthogonal  coordinate  system: 


D 

i! 

U 11  » 

0  g22 

0  \ 

0  r 

/  rnk. 

(g  )  = 

fi/gn 

0 

0 

1  'g2  2 

0  ' 
0 

^0  0 

g33/ 

1° 

0 

1/g33 j 

(B.  12) 


(B.  13) 


k 

using  these  in  (B.  13)  and  (B.  8)  we  can  write  the  matrix  of  d  ^  in  physical 
components  as: 


(d<k) ' 


/ 


1_ 
i  2 


;i 


\  — 
\2 
\ 


/g22  2  L  fu  1 

■ -  Y  -f  t -  Y  ^ 

<s7  \lg22  ’2 

/g33  3  /^TM 


gll  ;1  J  g33  ’3 


ell  1 

! -  \r 


ig  1 1  1  . 

—  v  „  + 
g22  ’2 


v 


;2 


/g33  1  3 
/ -  v  „  + 

\l  g22  ;2 


I  § 


22  2 
V 


gll  ^ 


/gl 1  1  ,  /g33  3 

\/  g33  ;3Jgll  ' 


I i 


lg  22  2 

-  V  0 

g33  ;3 


g22  2 
■—  v  „  + 

\l  g33  ’3 


;3 


IS 


33 


g 


22 


Using  (B.  6) 


(n(k))  =  (h1  \lgll  r2  ^g22  °3  \/g33 


Putting  these  last  two  in  (B.  11): 

d(n)  =  v:1l(nl>2«ll  +  v;2<n2)2  ®22  +  v;3<n3>2  g33  +  (g22  v ;1  +  «1 1  v;2) 


1  2 
n  n 


+ 


/  1  .  3\  13,/  3  ,  2\  2  3  rn  ia.\ 

Vgll  v;3+  %3  v;l)n  n  +  \*33  v;2  +  g22  v;3/  n  n  '  (B'14) 
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Since  we  desire  d^^  in  the  direction  of  acceleration,  we  will  express 
n  as  follows: 

_(\lg  11  A  \lg 22  A  \fg33  A*  \ 

lAj  Tat-  ) 

where 

|  A|  =  ygn(Ab2  +  g22(A^)2  +  g33(A3)2  =  acceleration, 
then  we  may  write  (B.14)  as: 


'(a) 


.  ^  <A'>2gn  +  4(a2>2  ^22  *  y:3(A3>2  «33  +  (S22  **  +  *11  4 >A‘  A' 


gjjtA1)2  +  g,,(A2)2  +  g,,<A3)2 


=22 


533 

2  4- 

V  ~  + 


22  v;3  g33 


V£) 


A2  A3 


g 


1 1 


(Al)2  +  S22(a2)2  +  ^^(a3)2 


(B.  15) 


’33 


For  a  rectangular  coordinate  system: 


1  2  3 

x  =  x,  x  =  y,  x 

gl  1  "  g22  g33  =  1 


z 


V  =  dx\/d  X1'  . 

;;i  ' 

For  a  cylindrical  coordinate  system: 


1 

X  -  r. 

2  -  n  3 

x  -  6,  x 

=  z 

gl  1  "  g33  1  ’  g22 

2 

-  r 

1  8r 

j  dr 

9 1-  ro 

ae  6 

dr 

3z 

i  j 

36  6 

36  r 

86 

v  .  -  j 

;j  ! 

dr  r 

36  r 

9z 

t 

1 

dz 

dz 

3z 

\  dr 

36 

3z 

For  a  spherical  coordinate  system: 

1  2  _  fl  3  _  . 

x  -  r,  x  -6,  x  -ip 

*11  *  *'  *22  *  r'2'  *33  *  r2  sin2<> 
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<V;J>  ■ 


/  M 

3r  ■ 

ss  - re 

3r 

/  9r 

3  ip 

|»  +  » 

|»  +£ 

30 

9r  r 

i 

90  r 

dip 

|f  +  ip  cot 0 

d± 

dip 

•  2, 


-  rip  sin  0 

-  ^sin0  cos0 


\ 


0cot0 


J 


Using  these  in  (B.  15)  we  can  write: 

For  a  rectangular  coordinate  system  in  three  dimensions: 


(A) 


£<a  )2 

ox  X 


+  |2(A  )‘ 


+ 


(A  )2 
3z  z 


+ 


3x  ,  3y 

a £+ a5 


A  A  +(|S  +  |z)a  a+(?2  +  |^  A  A 

x  y  \3z  3x  /  x  z  13z  9y;  y  z 


<AX)2  + 


(Ay)2  + 


<V 


For  a  cylindrical  coordinate  system: 


(A  )2  +  +  r  j  /  »  x2  +  9z  (A  x2  +  far  2  90 

9r  '  r  +  l  90  r  /  r  A0'  +  9z  z'  1  90  +  r 


J(A) 


a  ,  A  Afl 
3r  /  r  0 


(Ar)2  +  (rA0)2  +  (Az)2 


9r  jlz 
k9z  9r 


)  ArAz  +  ill  +  r2  H  i 


AflA 
0  z 


(Ar)2  +  (rAg)2  +  (Az)2 


For  a  spherical  coordinate  system: 


’(A) 


9r 


7  (A rf  + 


90  +  r 


:  )  0" A@)2  h  (jf  +  ^  ^  cot0)  (rA;  sin0)2 


(A^)2  +  (rA^)2  +  (rA^sin0)2 


iL 


(||+r2|^!  A  A.+ 


9r 


r  0 


9r,  2  .  2  Q  9i p  )  *  a  , 

a~  +  r  sin  0  A  A,  + 
dip  or  J  r 


+  l^+  2  Ml  r2  A  A 

‘ip  [dip  _  90]  A0<// 


(A^)2  +  (rA^)2  +  (rA^sin0)2 


It  is  easily  seen  from  the  above  that  for  a  one-dimensional  system: 


9x 

9x 


which  is  the  quantity  used  by  von  Neumann  to  activate  his  "q"  and  whose 
magnitude  was  used  is  forming  the  "q." 
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It  will  be  noted  that  any  direction  could  be  selected  for  n  and  in  case 

some  other  direction  is  desired,  (B.  14)  can  be  used  to  calculate  d,  .. 

(n) 

Quite  obviously  the  first  two  expressions  above  for  can  be  written 
as: 


9x  '*  *2  ■  ®i(A  )2  +  /®i  + 


\A) 


TT  <A  f  + 

ox  x  9y 


(Ax)2  +  (Ay)2 


lw  +_M 


A  A 

x  y 


if  cylindrical  symmetry  is  assumed  with  no  angular  motion  and  y  is  either  the 
radial  direction  or  the  other  rectangular  direction.  In  this  case  the  com¬ 
ponents  of  the  acceleration  vector  are  obtained  from  the  motion  equations  as 
given  in  Ref.  2,  p.  20  and  are  reproduced  below: 


fd  2 


A  =  -I 

x  p  I 


XX 


dT 


T 


i*\ 


3x 


9y  y  i  / 


,  /a t  as  s 
A  =!  *Z+  -JAL  +  ,.yy. 

y  p  I  ax  ay  y 


T. 


-  acceleration  in  x-direction 

n 

-  acceleration  in  y-direction 


ueo 


where 


2. .  r  normal  stress  in  i  direction 

r^(i  f  j)  £  shear  stress  in  the  j  direction  on  a  plane  with  a  normal 
in  the  direction  i 
p  =  mass  density . 


In  connection  with  the  HEMP  code  (see  Ref.  2)  the  expression  d^n)  has 
another  application.  This  is  in  connection  with  the  problem  involving  a 
material  with  a  different  yield  strength  in  compression  than  it  has  in  tension. 
In  such  a  case  one  must  be  able  to  determine  whether  there  is  compression  or 

tension  in  the  direction  of  motion.  The  quantity  d,  ,  can  answer  this  if  n  is 

In;  ^ 

along  the  velocity  vector.  ~ 


These  terms  are  set  to  zero  for  the  plane  case. 
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